function res = CampbellSchillerRegression_riskAdjust(Data,termPrem,fwdPrem,m,numBoots)

% The CS regression reads:
% y_t+m_(k-m) - y_t_k - term_t+m_(k-m) + term_t_(k-m) + m/(k-m)*fwdPrem_t_(k-m) ...
%   = alpha + betta*m/(k-m)*(y_t_(k)-y_t_(m)) + eps_t+1
%
% INPUT: Must be a data matrix of dim T * numYields, 
% where yields have maturities of 1 quarter, 2 quarters,3 quarters, ...
[T,n]          = size(Data);
matSelect      = m:m:n;
shorRate       = Data(1:T-m,m);
yieldsUsed     = Data(:,matSelect);
numYields      = size(yieldsUsed,2);
termPremUsed   = termPrem(:,matSelect);
fwdPremUsed    = fwdPrem(:,matSelect);
CSBetta        = NaN(2,numYields);
CSBetta_se     = NaN(2,numYields);
CSBetta_seBoot = NaN(2,numYields);
Ydata          = NaN(T-m,numYields);
Xdata          = NaN(T-m,numYields);
R2             = NaN(1,numYields);
for i=1:numYields-1
   Y  = yieldsUsed(1+m:T,i)-yieldsUsed(1:T-m,i+1) ...
       - termPremUsed(1+m:T,i)...
       + termPremUsed(1:T-m,i)...
       + fwdPremUsed(1:T-m,i)*m/(matSelect(1,i+1)-m);
   %x  = (yieldsUsed(1:T-m,i+1)-shorRate)*m/(matSelect(1,i+1)-m);
   x  = yieldsUsed(1:T-m,i+1)-shorRate;
   X  = [ones(T-m,1) x*m/(matSelect(1,i+1)-m)];
   resOLS            = nwest(Y,X,m+1);
   CSBetta(:,i+1)    = resOLS.beta;
   CSBetta_se(:,i+1) = resOLS.beta./resOLS.tstat;
   Ydata(:,i+1)      = Y;
   Xdata(:,i+1)      = x*m/(matSelect(1,i+1)-m);
   R2(:,i+1)         = resOLS.rsqr;
   if numBoots > 0
       % Bootstrapping the SE
       bootWidth = 10;
       %numBoots  = 1D5;
       rng(1,'twister')
       bsdata=block_bootstrap([Y,X],bootWidth,numBoots);
       betaBoot  = zeros(2,numBoots);
       for s=1:numBoots
           dataBoot = bsdata(:,:,s);
           Yboot    = dataBoot(:,1);
           Xboot    = dataBoot(:,2:3);
           betaBoot(:,s) = (Xboot'*Xboot)\(Xboot'*Yboot);
       end
       CSBetta_seBoot(:,i+1) = std(betaBoot,0,2);
   end
end

% Output
res.CSbetta        = CSBetta;
res.maturities     = matSelect;
res.CSbetta_se     = CSBetta_se;
res.CSBetta_CI95   = [CSBetta(2,:)-1.96*CSBetta_se(2,:);CSBetta(2,:)+1.96*CSBetta_se(2,:)];
res.CSBetta_CI99   = [CSBetta(2,:)-2.575*CSBetta_se(2,:);CSBetta(2,:)+2.575*CSBetta_se(2,:)];
if numBoots > 0
    res.CSBetta_seBoot = CSBetta_seBoot;
    res.CSBetta_CI95Boot =[CSBetta(2,:)-1.96*CSBetta_seBoot(2,:);CSBetta(2,:)+1.96*CSBetta_seBoot(2,:)];
end
res.Ydata          = Ydata;
res.Xdata          = Xdata;
res.R2             = R2;


